#ifndef _COMMON_H_

#define _COMMON_H_

#ifdef BOOT_FOUND
  #include <boost/numeric/ublas/matrix.hpp>
  #include <boost/numeric/ublas/vector.hpp>
  #include <boost/numeric/ublas/matrix_proxy.hpp>
  #include <boost/numeric/ublas/io.hpp>
  using namespace boost::numeric::ublas;
#else
  #include<tvmet/Matrix.h>
  #include<tvmet/Vector.h>
  /// This is just to ensure that I can switch b/w boost and tvmet
  /// to measure performance
  #define DECL_MAT( type , dim ) \
          tvmet::Matrix<type,dim,dim>
  #define DECL_VEC( type , dim ) \
          tvmet::Vector<type,dim>
#endif

// Some standard libs
#include<vector>
#include<fstream>
#include<iostream>
#include<blitz/array.h>
#include<blitz/tinyvec.h>

#define DECL_UD_VEC( name , T , dim ) \
template< class T , int dim > \
  struct name { \
    T data[dim]; \
    inline name<T,dim> operator+( name<T,dim> &other ) { \
      name<T,dim> result = *this ; \
      for( int i = 0 ; i < dim ; i++ ) \
        reinterpret_cast< T * >(&result)[i] += reinterpret_cast< T * >( &other )[0] ; \
      return result; \
  } \
  inline name<T,dim> operator-( name<T,dim> &other ) { \
    name<T,dim> result = *this ; \
    for( int i = 0 ; i < dim ; i++ ) \
      reinterpret_cast< T * >(&result)[i] -= reinterpret_cast< T * >( &other )[0] ; \
      return result; \
  } \
  inline name<T,dim> operator*( T other ) { \
    name<T,dim> result; \
    for( int i = 0 ; i < dim ; i++ ) \
      reinterpret_cast< T * >(&result)[i] *= other ; \
      return result; \
  } \
};\
template<class type , int dim>\
  std::ostream& operator<<(std::ostream &out , const name<type,dim> &res) {\
  out << "[ ";\
  for( int i = 0 ; i < dim ; ++i )\
    out << res.data[i] << " "; \
  out << "]";\
  return out;\
}


#define MULTICOM_BZ_UD_VEC( name , type , dim ) \
namespace blitz { \
  template<> struct multicomponent_traits<name<type,dim> > { \
    typedef type T_element; \
    static const int numComponents = dim; \
  }; \
}

DECL_UD_VEC( point , dtype , dim )

MULTICOM_BZ_UD_VEC( point , double , 2 )
MULTICOM_BZ_UD_VEC( point , double , 3 )
MULTICOM_BZ_UD_VEC( point , double , 4 )
MULTICOM_BZ_UD_VEC( point , double , 5 )
MULTICOM_BZ_UD_VEC( point , double , 6 )

MULTICOM_BZ_UD_VEC( point , float , 2 )
MULTICOM_BZ_UD_VEC( point , float , 3 )
MULTICOM_BZ_UD_VEC( point , float , 4 )
MULTICOM_BZ_UD_VEC( point , float , 5 )
MULTICOM_BZ_UD_VEC( point , float , 6 )

#define DECL_BZ_VEC( type , dim ) \
        blitz::Array< point<type,dim> , dim >
#define DECL_BZ_VEC_N( type , dim , N ) \
        blitz::Array< point<type,N> , dim >
#define DECL_BZ_VEC_NxN( type , dim , N ) \
        blitz::Array< point<type,N,N>, dim >
#define DECL_TVEC( type , dim ) \
        blitz::TinyVector<type,dim>
#define DECL_BLOCK_VEC( itype , dtype , dim ) \
        std::vector< pavan::block<itype,dtype,dim> >
#define DECL_BLOCK( itype , dtype , dim ) \
        pavan::block<itype,dtype,dim>
#define DECL_ZONE_VEC( itype , dtype , dim ) \
        std::vector< pavan::zoneInterface<itype,dtype,dim> >
#define DECL_ZONE( itype , dtype , dim ) \
        pavan::zoneInterface<itype,dtype,dim>

#define DECL_SCALAR_DATA( dtype , dim ) \
        std::vector< blitz::Array<dtype,dim> >

#define DECL_VECTOR_DATA( dtype , dim , N ) \
        std::vector< DECL_BZ_VEC_N( dtype , dim , N ) >

#define DECL_LOOP_2( var , itype , dim )  \
    for( int ivar = 0 ; ivar < dim ; ++ivar ) \
      for( itype i = var.lbound(0) ; i <= var.ubound(0) ; ++i ) \
        for( itype j = var.lbound(1) ; j <= var.ubound(1) ; ++j )

#define DECL_LOOP_3( var , itype , dim )  \
    for( int ivar = 0 ; ivar < dim ; ++ivar ) \
      for( itype i = var.lbound(0) ; i <= var.ubound(0) ; ++i ) \
        for( itype j = var.lbound(1) ; j <= var.ubound(1) ; ++j ) \
          for( itype k = var.lbound(2) ; k <= var.ubound(2) ; ++k )

struct Time {
  Time() : _startTime(0.0) , 
           _endTime(0.0) ,
           _deltaT(0.0) ,
           _iter(0) , _curIter(0) { }

  double _startTime;
  double _endTime;
  double _deltaT;
  long unsigned int _iter;
  long unsigned int _curIter;
};

static const char _nameIJK[] = { 'I' , 'J' , 'K' };

#define DEFAULT_PADDING 2
/*namespace blitz {

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::transposeSelf( TinyVector<int,N_rank> &perm  )
{
    int dummy = 0;
    for( int i = 0 ; i < N_rank ; ++i ) dummy += perm[i];
    BZPRECHECK( dummy == N_rank * (N_rank-1) / 2 ,
        "Invalid array transpose() arguments." << endl
        << "Arguments must be a permutation of the numerals (0,...,"
        << (N_rank - 1) << ")");

    // Create a temporary reference copy of this array
    Array<T_numtype, N_rank> x(*this);

    // Now reorder the dimensions using the supplied permutation
    for( int i = 0 ; i < N_rank ; ++i )
      doTranspose( i , perm[i] , x );
}


}*/

#endif
